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IMPLEMENTATION  OF  A  BALANCE  OPERATOR  IN  NCOM 


1.  INTRODUCTION 
1.1  Balance  Constraints 

Due  to  the  sparsity  of  oceanographic  observations,  a  common  practice  in  oceanic  data  assimila¬ 
tion  is  to  constrain  the  model  state  increments  to  be  in  (approximate)  geostrophic  and  hydrostatic 
balance.  This  strategy  allows  one  to  extract  the  most  persistent  modes  of  variability  from  the 
data  and  filter  out  smaller /faster  scale  motions  that  are  barely  resolved  by  the  observations.  The 
geostrophic  and  hydrostatic  balance  constraints  are  supplemented  by  the  (linearized)  equation  of 
state  and,  quite  often,  by  additional,  linear  relationships  between  temperature  and  salinity  (T  —  S 
constraints),  and/or  an  integral  continuity  constraint. 

The  basics  of  the  above  mentioned  “balance  operator  technique”  were  originally  developed  for 
meteorological  data  assimilation  (Derber  and  Bouttier  1999;  Cullen  2003)  and  later  implemented 
in  oceanographic  applications  of  variational  data  assimilation  (Weaver  et  al.  2005). 

To  ensure  computational  efficiency,  the  balance  constraints  are  introduced  as  a  sequence  of 
linear  operations  on  the  “unbalanced”  (statistically  independent)  components  of  the  ocean  state 
xi  =  [T(x,  z),  S(x,  z)],  which  produces  the  balanced  constituents  of  the  remaining  part  of  the  state 
vector  xo  =  [u(x,  z),  £(*)].  After  that,  the  unbalanced  constituents  X2  of  X2  are  added  to  obtain  the 
full  update  of  the  ocean  state  vector.  To  be  more  specific,  assume  that  the  model  state  variables 
{T,  S',  C,  u}  are  represented  by  the  M-dimensional  vector  y  €  1ZM  of  the  model  fields’  values  at 
the  computational  grid  points.  The  balance  constraints  are  introduced  by  partitioning  y  as  follows 
(Weaver  et  al.  2005): 


y  = 


Xl 

X2 


Xl  = 


e  n 


Mi 


X2  = 


e  n 


m2 


Mi  +  M2  =  M, 


(1) 


and  representing  X2  in  the  form  X2  =  Lxi  +  X2,  where  L  stands  for  the  balance  operator  and 
X2  denotes  the  unbalanced  constituents  of  X2.  The  discretized  version  of  the  operator  L  can  be 
symbolically  represented  by  an  Mi  x  M2  matrix 


eq.of  state 

u 

geostrophy 

T 

_c_ 

hydrostatics 

S 

continuity 

which  is  a  finite-difference  discretization  of  the  following  constraints  (in  the  sequence  they  are 
symbolically  written  above): 
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p  =  P0  +  a(x,  z)T  +  P(x,z)S , 

(3) 

u=  ft  \  k  x 

f(x)po 

(4) 

dzp  =  - gp ,  p(x,  0)  =  p0g((x), 

0 

(5) 

0  =  div  J  udz, 

(6) 

h(x) 


where  x  =  (x,  y )  is  the  horizontal  coordinate,  z  is  the  vertical  coordinate,  po  is  the  mean  density  of 
seawater,  a  and  (3  are  the  Taylor  expansion  coefficients  in  the  linearized  equation  of  state  (which 
depend  on  the  spatial  coordinates  because  the  background  temperature  Tj,  and  salinity  fields 
do),  /  is  the  Coriolis  parameter,  k  is  the  vertical  unit  vector,  V  is  the  horizontal  gradient,  p  is 
the  pressure,  h  is  the  ocean  depth,  and  g  is  the  gravitational  acceleration.  In  applications,  these 
constraints  are  usually  modified  to  form  the  following  sequence  of  operations  which  linearly  map 
the  input  xi  =  [T,  S]  to  the  output  X2  =  [u,  (]: 


P  =  Po  +  «(*,  z)T  +  /3(x,z)S , 


0  0 


Vh(x)V C  =  div  J  j  V p(x,  z')dz'dz, 

h(X)  * 

0 

PoC(x)  +  J  p(x,  z')dz' 


P  =  9 


u  = 


f(x)po 


k  x  Vp. 


(7) 

(8) 

(9) 

(10) 


The  equation  (8)  involves  an  iterative  inversion  of  the  elliptic  operator  V/iV,  which  increases  the 
computational  cost  of  L.  For  that  reason,  the  respective  constraint  (integral  continuity)  is  sometimes 
replaced  by  assigning  “a  level  of  no  motion”  (setting  p  =  0)  at  a  certain  reference  depth  zref,  or  at 
the  bottom: 

0 

C  =  -—  f  p(x,z)dz.  (11) 

Po  J 

h 

A  more  general  version  of  L  is  obtained  by  splitting  the  salinity  field  perturbations  S  into  balanced 
S  and  unbalanced  S  components: 

S  ->S  +  S  =  6(x,z)T  +  S,  (12) 

where  6  is  the  user-supplied  scalar  field  describing  the  spatial  variability  of  the  T  —  S  relationship. 


Since  the  balanced  and  unbalanced  components  are  assumed  to  be  uncorrelated,  it  is  desirable 
to  preserve  the  total  salinity  error  variance  estimated  from  the  previous  analyses.  In  the  presented 
formulation  of  the  balance  operator,  this  constraint  is  introduced  by  the  user-defined  coefficient 
7  (x,z): 

S  =  ^S+  ^1  -y2  S  =  'yOT  +  ^1  -72  S.  (13) 
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In  the  upper  layers  of  the  ocean,  where  temperature  and  salinity  perturbations  are  weakly  corre¬ 
lated,  7  is  usually  set  to  zero,  while  in  the  deep  layers  it  is  useful  to  set  it  fairly  close  to  1  (Ricci 
et  al.  2005).  From  the  viewpoint  of  numerics,  introduction  of  the  balanced  salinity  component  is 
reduced  to  redefinition  of  the  coefficients  a  and  (3  in  the  equation  of  state: 

a — >  a( l  +  'yO),  (3 — ►  (3\j  1  -  y2.  (14) 


In  the  present  report,  we  describe  the  numerical  implementation  of  Eq.  (7)— (12) ,  which  has 
been  derived  from  the  NCOM  code  (Section  2.1).  Since  the  adjoint  of  L  is  crucial  for  constraining 
the  background  covariance  (BEC)  and  the  descent  process  used  in  the  data  assimilation,  the  code 
for  Lt  has  been  also  constructed  (Section  2.2). 


1.2  Adjoint  of  the  Balance  Constraints  and  the  Background  Error  Covariance 


Partitioning  (1)  of  the  state  vector  implies  the  following  structure  for  the  BEC  matrix  (assuming 

(xiX2)  =  0): 


(xlxl~) 

(xi4) 

ER  BlLT 

E!  0 

Bi  0 

Ex  Lt 

_(X2XI) 

(X2X2)_ 

LBi  LBiLt  +  B2 

L  E2 

0  b2 

0  e2 

Here  Bi  =  (xixi)  and  B2  =  (x2x2)  are  the  M \  x  M\  and  M2  x  M2  BEC  matrices  of  the  unbalanced 
components  of  y  and  E]  2  are  the  identity  matrices  of  the  respective  sizes. 


In  many  applications  (such  as  NCODA  3dVar  and  4dVar),  Bij2  are  represented  by 

Bi,2  =  Vi)2Ci,2Vi>2,  (16) 

where  V i>2  are  the  diagonal  rms  background  error  variance  matrices  of  xi  and  x2,  and  Cg2  are 
the  respective  correlation  matrices  modeled  by  the  polynomials  of  the  diffusion  operator  A  (e.g., 
Weaver  and  Courtier  2001;  Yaremchuk  and  Sentchev  2012;  Yaremchuk  et  al.  2013).  In  particular, 
the  correlation  models  currently  available  at  NRL  have  the  form: 

Cp-  [E-  ■jAp2,  Ce~exp[^-A],  (17) 

where  a  and  b  are  the  decorrelation  radii  and  A  is  the  Laplacian  operator.  Note  that  Cp  can  be 
used  in  both  state-space  and  observation  space  3d/4dVar  formulations,  whereas  Cp  is  not  invertible 
and,  therefore,  can  be  used  either  in  an  observation  space  formulation  or  in  a  B-preconditioned 
state-space  approach. 


Substitution  of  (16)  into  (15)  provides  the  following  factorizations  of  B  and  B 


-l. 


where 


B-1  =  V-1C-1V-T,  B  =  VCVt, 


o 

u 

1 _ 

i 

< 

o 

[  L  E2 

- 1 

CM 

> 

o 

_ l 

C 


C,  0 

0  c2 


(18) 


(19) 
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Furthermore,  since 


1  _ 

fvr1 

1 

o 

i - 

o 

UJ 

1 _ 

o 

1  CM 

> 

1 

1 

r~ 

m 

to 

l _ 

(20) 


the  balanced  BEC  and  its  inverse  allow  explicit  symmetric  factorization  using  the  square  roots  of 
the  correlation  matrices  (17): 

B  1  =  B~t/2B^1/2,  B  =  B1/2Bt/2,  with  B1/2  =  VC1/2,  (21) 


which  is  an  important  property  ensuring  the  computational  efficiency  of  many  data  assimilation 
algorithms.  Expressions  for  the  balanced  BEC  and  its  inverse  have  the  form: 


i  _ 

Bj-1  +  LtB71L 

-LTB2r 

R  _ 

B!Lt 

B2  1 

5  K  — 

LBi  LBr1LT  +  B2 

(22) 


which  explicitly  shows  the  necessity  of  having  the  adjoint  code  (LT)  for  operations  with  the  balanced 
BEC  and  its  inverse. 

Finally,  the  matrix  G  =  WT  can  be  effectively  used  as  a  natural  metric  in  the  space  of 
cost  function  gradients.  The  associated  geometry  inhibits  descent  in  the  unbalanced  directions 
and,  therefore,  has  the  potential  to  make  the  3d/4dVar  optimization  processes  more  efficient  in 
recovering  balanced  corrections  to  the  background  state. 

2.  CODE  DESCRIPTION 

The  algorithms  that  comprise  the  NCOM  computational  kernel  are  described  in  detail  by 
Martin  (2000),  with  some  of  the  more  recent  additions  described  by  Morey  et  al.  (2003)  and  Barron 
et  al.  (2006).  The  vertical  mixing  models  that  are  used  in  NCOM  are  the  Mellor-Yamada  level  2 
(Mellor  and  Yamada  1974)  and  level  2^  (Mellor  and  Yarnada  1982)  turbulence  closure  schemes.  The 
equation  of  state  used  is  that  of  Mellor  (1991).  There  are  options  for  several  advection  schemes, 
of  which  the  most  commonly  used  is  third-order  upwind  (Holland  et  al.  1998),  which  implicitly 
includes  biharmonic  diffusion. 

The  balance  operator  is  written  as  a  subroutine,  ncom_balanced,  which  is  called  to  compute 
the  sea-surface  height  (SSH)  and  velocity  perturbations  from  the  input  temperature  and  salinity 
perturbations.  As  noted  in  the  previous  section,  two  methods  of  computing  the  SSH  perturbations 
are  provided:  (a)  the  solution  of  a  elliptic  equation  formed  from  the  horizontal  divergence  of  the 
linearized,  depth-integrated,  barotropic  and  baroclinic  pressure  gradient  terms  of  the  horizontal 
momentum  equations  and  the  depth-integrated  continuity  equation  (8),  and  (b)  a  dynamic  height 
calculation  using  a  specified  level  of  no  motion  (Eq.  (11)). 

2.1  Dynamic  Height  Calculation  for  SSH  Perturbations 

Subroutine  denp_lin  is  called  to  compute  density  perturbations  p'  from  temperature  and  salin¬ 
ity  perturbations,  T'  and  S' ,  respectively,  using  the  linearized  equation  of  state 

p  =  a(x,  z)T'  +  /3(x,  z)S' .  (23) 
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The  equation  of  state  is  linearized  around  reference  values  of  temperature  and  salinity  using  coef¬ 
ficients  of  thermal  and  salinity  expansion,  a  and  (3,  respectively. 


Subroutine  bpg_term  is  called  to  compute  baroclinic  pressure  gradient  perturbation  terms  from 
the  density  perturbations  p'.  These  are  computed  using  the  same  numerics  as  used  in  NCOM 
(Martin  2000),  i.e. ,  on  the  sigma  coordinate  part  of  the  grid,  the  baroclinic  pressure  gradient 
perturbation  in  the  x-direction  is  calculated  for  the  kth  layer  as 


9 


1  dp' .  1  dp' . 

- —\k  = - —  L_i  + 

p0  dx  p0  dx  p0 Ax1 

1 


(-L)“(Ac7fc_i  +  A  ak)5x(p'k_1  +  p'k ) 


-  -(o-fc-i  +  crk)(5xD)(Szpr  )), 


(24) 


and  on  the  z-level  part  of  the  grid,  the  baroclinic  pressure  gradient  perturbation  is  calculated  as 

1  dp'  1  dp'  g  1 

p0  dx  k  p0  dx  k  1  p0Axv 


r [,Azk—\8xpk_ l  T  Azk5xpk), 


(25) 


where  Axu  is  the  local  grid  spacing  in  the  x  direction  at  a  u  point,  Aak  is  the  fractional  thickness 
of  sigma  layer  A;  at  a  T  point,  A zk  is  the  layer  thickness  on  the  z-level  part  of  the  grid,  Sx  is  the 
differential  operator  in  the  x  direction,  Sz  is  the  differential  operator  in  the  vertical  direction,  Du 
is  the  total  depth  of  the  sigma  part  of  the  grid  at  a  it  point,  D  is  the  total  depth  of  the  sigma  part 
of  the  grid  at  a  T  point,  and  the  overbar  indicates  a  horizontal  average  in  the  specified  direction. 
The  baroclinic  pressure  gradient  perturbation  in  the  y-direction  is  calculated  analogously. 


Subroutine  dynJit  is  called  to  compute  SSH  perturbations  from  the  density  perturbations 
p'  using  a  dynamic  height  calculation  and  an  assumed  level  of  no  motion,  i.e., 

C'  =  --5>fcA*fc,  (26) 

Po  ^ 

where  the  sum  is  computed  from  the  assumed  level  of  no  motion  to  the  surface.  Note  that  the  layer 
thickness  for  layers  on  the  sigma  part  of  the  grid  is  computed  as  A  zk  =  DAak. 


Subroutine  vel_geo  is  called  to  compute  velocity  perturbations  u',  v'  from  the  baroclinic  pres¬ 
sure  gradient  and  SSH  perturbations  assuming  geostrophy,  i.e., 


J_J _dj/_  _  g  5y('xy 

Jy  Po  dy  Jy  A yv 


(27) 


_1  1  dp'  g  Sx  CXV 

Jx  po  dx  Jx  Axu 


(28) 


where  /  is  the  local  Coriolis  parameter,  A yv  is  the  local  grid  spacing  in  the  y  direction  at  a  v  point, 
and  the  overbars  indicates  horizontal  averaging  in  the  specified  direction(s). 


2.2  Solution  of  an  Elliptic  Equation  for  the  SSH  Perturbations 


The  alternate  method  for  computing  the  SSH  perturbations  is  to  solve  Eq.  (8),  an  elliptic 
equation  for  the  SSH  perturbations.  The  finite  difference  form  of  this  equation  that  is  used  (which 
is  similar  to  the  form  used  within  NCOM  to  update  the  SSH)  is 
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Fig.  1  Flow  chart  for  the  balance  operator  code  (left),  and  for  the  adjoint  code  for  the  balance  operator  (right). 


4( 


AyugHu 

Axu 


SxC)  +Sy( 


A  xvgHv 
A  yv 


Sy?) 


-5x(AyuJ2 


A  z%  dp' 

Po  dx 


k) 


A  z\  dp' 

Po  dy 


k) 


(29) 


where  Hu  is  the  total  depth  at  a  u  point,  Hv  is  the  total  depth  at  a  v  point,  and  the  sums 
are  taken  in  the  vertical  over  the  whole  depth  of  the  water  column.  This  equation  is  solved  for 
the  SSH  perturbations  using  a  pre-conditioned,  conjugate-gradient  method.  Note  that,  relative  to 
Eq.  (8),  in  forming  Eq.  (29),  the  x-  and  y  component  equations  were  multiplied  by  the  grid-cell  area 
before  taking  the  horizontal  divergence  so  that  the  coefficients  input  to  the  elliptic  solver  would  be 
symmetric. 


2.3  Direct  Code 


The  flow  chart  on  the  left  side  of  Figure  1  shows  the  structure  of  the  code  for  the  balance 
operator  realized  in  subroutine  ncom_balanced.  The  NCOM  grid  parameters  are  passed  into 
ncom_balanced  in  the  subroutine  argument  list  to  improve  the  flexibility  of  the  code  in  accom¬ 
modating  nested  grids.  The  subroutine  sequentially  solves  the  system  of  Eq.  (T)-(ll)  discretized  in 
a  manner  consistent  with  the  NCOM  formulation. 

Specifically,  the  constituents  of  ncom_balanced  perform  the  following  operations: 

•  denp_lin:  computes  density  perturbations  via  Eq.  (7), (23) 

calls  ce_mel3  to  compute  expansion  coefficients  a  and  (3  using  7  and  6  via  Eq.  (14). 

•  bpg_term  -  computes  baroclinic  pressure  gradients  via  Eq.  (9),(24)-(25). 
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•  cgssor:  solves  for  SSH  using  Eq.  (8), (29) 

calls  minmaxl 

•  dynJit:  computes  SSH  using  Eq.  (11), (26) 

•  vel_geo:  computes  horizontal  velocities  using  Eq.  (10),(27)-(28) 

calls  avev  and  aveu  to  interpolate  between  the  u  and  u-grids. 

Switching  between  the  modes  is  made  by  changing  the  parameter  iab,  i.e. ,  iab=l  invokes  the 
b-mode  and  iab=2  invokes  the  a-mode. 

2.4  Adjoint  Code 

The  flow  chart  on  the  right  side  of  Figure  1  shows  the  structure  of  the  code  for  the  adjoint 
of  the  balance  operator  realized  in  subroutine  ncom_balanced.  Similar  to  ncom_balanced,  the 
grid  parameters  are  passed  into  ncom_balancea  within  the  subroutine  argument  list  to  improve 
the  flexibility  of  the  code  in  accommodating  nested  grids.  The  subroutine  sequentially  solves  the 
adjoints  of  Eq.  (T)-(ll)  in  the  reverse  order. 

Specifically,  the  components  of  ncom_balancea  perform  the  following  operations: 

•  vel_geoc:  computes  the  adjoint  of  the  velocity  perturbation  Eq.  (10),(27)-(28) 

calls  avevc,  aveuc,  the  adjoints  of  the  interpolation  between  the  u  and  v  grids. 

•  cgssorc:  solves  the  adjoint  of  Eq.  (8), (29) 

calls  minmaxl. 

•  dynJit  c:  computes  the  contribution  to  the  adjoint  density  perturbation  Eq.  (11),(26). 

•  bpg_termc:  computes  the  contribution  to  the  adjoint  density  from  the  adjoint  pressure  gradients 
Eq.  (9),(24)-(25). 

•  denp_linc:  computes  the  adjoint  temperature  and  salinity  Eq.  (7), (23). 

calls  ce_mel3  to  compute  expansion  coefficients  a  and  (3  using  7  and  9  Eq.  (14). 

3.  TESTING 

3.1  NCOM  Configuration 

The  model  was  configured  at  3-krn  resolution  on  an  85x294  horizontal  grid,  with  32  levels  in  the 
vertical.  The  top  22  a  levels  follow  the  bathymetry  from  the  surface  to  a  maximum  depth  of  291  nr, 
and  10  fixed-depth  levels  are  used  below  291  nr.  Initial  and  open  boundary  conditions  for  the  SSH  (, 
temperature  T,  salinity  S,  and  horizontal  velocities  u,v  were  obtained  from  Global  NCOM,  which 
was  run  operationally  at  the  Naval  Oceanographic  Office  (Barron  et  al.  2004).  Tidal  forcing  was  not 
used  in  this  application.  The  Adriatic  model  was  forced  by  river  runoff  and  by  atmospheric  fields 
derived  from  the  8-kin  horizontal  resolution,  regional  ALADIN  atmospheric  model  (Ivatek-Sahdan 
and  Tudor  2004). 

In  the  described  numerical  experiments,  the  state  vector  y  comprised  all  the  grid  point  val¬ 
ues  of  £,  T,  S,  u,  and  v.  With  the  given  3-dimensional  grid  and  bathymetry,  the  state  vector  has 
M=l,493,570  components.  The  dimension  M2  of  the  balanced  constituent  was  equal  to  743,526 
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20  40  60  80  20  40  60  80  20  40  60  80 

Fig.  2  —  Maps  of  the  diagonal  elements  of  V2  corresponding  to  the  unbalanced  rms  error  variance  of  the  horizontal 

velocity  at  three  er-levels. 

(the  number  of  u,C  grid  points).  The  background  values  of  y (t)  were  taken  from  the  NCOM  sim¬ 
ulation  described  by  Martin  et  al.  (2009)  and  then  adjusted  to  suppress  temperature  and  salinity 
biases  during  the  assimilation  time  interval  (from  0.00  UTC  on  08/14  to  0.00  UTC  on  08/29/2006). 
After  the  adjustment,  the  horizontally-  and  temporally-averaged  misfits  between  the  background 
solution  and  the  T  and  S  observations  did  not  exceed  0.02°C  and  0.005  psu,  respectively. 

3.2  Balanced  and  Unbalanced  Perturbations 


The  relative  magnitude  of  the  ageostrophic  motions  in  the  model  solution  was  assessed  by  the 
non-dimensional  parameter 


Kz(t)  =  R 


u.  /  -* \  | xyz 
|div(tt)| 

.  .  , xyz  5 

\ub\ 


(30) 


where  u(x,z,t)  is  the  horizontal  velocity  of  the  solution,  b  denotes  the  background  or  reference 
solution,  the  overbar  indicates  3d  averaging  in  the  depth  interval  between  the  surface  and  depth 
z,  and  R  =  9  km  is  the  Rossby  radius  of  deformation  (Cushman- Rosin  and  Korotenko  2007).  The 
definition  (30)  does  not  account  for  the  divergence  of  the  geostrophic  currents,  which  is  negligible 
in  the  considered  regional  application. 


To  assess  the  impact  of  the  increments  generated  by  the  data  assimilation  procedure,  we  per¬ 
turbed  the  initial  conditions  of  the  background  solution.  The  3d  structure  for  the  simulated  incre¬ 
ments  were  generated  as  follows.  First,  we  defined  the  background  error  correlation  model  Cp  using 
Eq.  (17).  Second,  300  eigenvectors  corresponding  to  the  300  largest  eigenvalues  of  B  were  computed, 
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Fig.  3  -  Maps  of  SSH,  temperature  (color)  and  velocity  perturbations  on  three  cr-levels  constrained  by  the  La  operator 


solving  the  generalized  eigenvalue  problem  B  1ei  =  A/G  1ei  and  the  last  (/  =300)  eigenvector  e 
was  elected  to  compute  the  perturbations. 

The  actual  perturbations  5x  of  the  initial  conditions  were  defined  by  rescaling  the  eigenvector 
by  a  factor  e  whose  magnitude  was  chosen  in  such  a  way  that  the  cr-grid  root-mean-square  value 
of  the  temperature  perturbation  field  was  equal  to  0.8°C. 

Sx  =  ee30o-  (31) 

The  diagonal  elements  of  Vi  and  V2  in  Eq.  (19)  were  obtained  as  time-averaged  rms  variances  of  the 
background  fields.  Since  the  diagonal  of  V2  represents  the  ageostrophic  (unbalanced)  constituents 
of  the  velocity  and  SSH  fields,  the  respective  background  rms  variances  were  multiplied  by  the 
mean  value  of  k  (Section  3.3),  which  characterizes  the  level  of  the  ageostrophic  motions  in  the 
background  solution. 

Three  types  of  perturbations  were  tested: 

•  balanced  perturbations  (Fig.  3)  computed  with  the  integral  continuity  constraint  (operator  La 
in  Eq.  (19)  defined  by  Eq.  (7)— (10)); 

•  balanced  perturbations  (Fig.  4)  computed  with  the  assumption  of  zero  pressure  anomalies  at 
the  bottom  (operator  L&  defined  by  Eq.  (7),  (8)— (11)) 

•  unbalanced  perturbations  computed  with  L  =  0  in  Eq.  (19). 
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Fig.  4  -  Maps  of  SSH,  temperature  (color)  and  velocity  perturbations  on  three  cr-levels  constrained  by  the  Lb  operator. 


3.3  Results 


The  background  run  of  the  model  was  executed  from  August  9  to  August  24,  2006.  Since  the 
Adriatic  Sea  is  mostly  shallow  (depths  less  than  280  m  occupy  78%  of  the  area),  the  value  of  k  was 
computed  by  averaging  over  the  upper  22  a  layers.  The  mean  value  of  k  during  the  background 
integration  period  was  0.38,  indicating  the  presence  of  a  significant  ageostrophic  component,  which 
was  primarily  associated  with  inertial  oscillations  and  upper-layer  Ekman  dynamics. 


The  time  evolution  of  the  difference  between  the  values  for  k  for  the  perturbed  and  background 
solutions  is  shown  in  Fig.  5.  As  can  be  seen  in  Fig.  5,  k  grows  by  a  factor  of  1.6  (from  0.4  to 
0.64)  after  the  perturbation  of  the  model  on  August  9  with  unbalanced  perturbations,  whereas  k 
experiences  only  modest  growth  for  the  6-balanced  (32%)  and  a-balanced  (19%)  perturbations. 


The  advantage  of  lower  ageostrophic  shocks  provided  by  the  a-type  perturbations  comes  at  an 
additional  computational  cost:  application  of  the  a-type  balance  operator  La  consumes  5.5  times 
more  CPU  time  than  application  of  L;,  (0.171  sec  against  0.031  sec  on  a  single  2.3  Ghz  processor). 
For  the  tested  state  vector  (M  =1,493,570),  these  CPU  requirements  can,  however,  be  considered 
to  be  negligible  compared  to  the  CPU  time  required  by  the  model  integration  for  1  day  (60  sec  on 
single  processor).  For  that  reason,  the  balance  operator  and  its  adjoint  can  be  considered  to  be 
relatively  computationally  inexpensive  components  of  an  assimilation  cycle. 
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time,  days 

Fig.  5  —  Time  evolution  of  the  difference  in  the  values  of  k  between  the  background  solution  and  the  three  types  of 

perturbed  solutions. 


4.  SUMMARY 

The  balance  operator  code  that  has  been  developed  can  be  used  in  NRL  data  assimilation 
systems  for  the  purpose  of  constraining  the  innovations  to  a  “slow”  (hydrostatically  and  geostroph- 
ically  balanced)  manifold.  This  will  allow  efficient  extraction  of  the  most  persistent  modes  from 
the  incoming  data.  More  importantly,  the  balance  operator  approach  will  provide  more  flexibility 
in  constraining  the  background  covariances  to  the  user-defined  subspaces  fully  consistent  with  the 
numerics  of  the  parent  model.  This  will  enable  a  flexible,  modular  approach  to  the  construction  of 
the  covariance  operators  and  require  code  development  for  the  dynamical  constraints  defined  by  the 
user  together  with  their  adjoints  (a  relatively  inexpensive  procedure  compared  to  the  modification 
of  the  existing  NCODA  system,  where  the  balance  constraints  are  hard-wired  in  the  correlation 
operator  code).  The  balance  operator  could  be  used  in  upgrading  NCODA  3d-  and  4dVar  (Ngodock 
and  Carrier  2014)  data  assimilation  systems  and  could  be  also  useful  in  similar  atmospheric  appli¬ 
cations  (Xu  and  Rosmond  2004;  Rosmond  and  Xu  2006). 

The  balance  operator  code  presented  in  this  report  can  be  further  developed  to  include  the 
parameterization  of  equatorial  dynamics,  non-linear  terms  on  the  right-hand  side  of  the  elliptic 
equation  for  the  SSH  perturbations,  and  MPI  parallelization. 
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